home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / MFLOAT10.ZIP / MFLOAT.DOC < prev    next >
Text File  |  1993-04-28  |  71KB  |  1,679 lines
  1.                MFLOAT Version 1.00
  2.                ===================
  3.  
  4. A brief description of "MFLOAT":
  5.  
  6. * "MFLOAT" handles arithmetic using high precision numbers of up to 72 decimal
  7.   digits.
  8.   
  9. * "MFLOAT" is optimized for maximum speed for high precision numbers.
  10.   These subroutines are the fastest for the 8086 microprocessor known to us.
  11.   If you find faster subroutines for a 8086, please send me a mail.
  12.  
  13. * "MFLOAT" subroutines are much faster than a standard emulation of a 
  14.   coprocessor for the same accuracy. "MFLOAT" does not use the coprocessor.
  15.  
  16. * "MFLOAT" is written in assembler for the 8086 microprocessor and
  17.   compatibles. The calculation is 10 to 20 times faster than that of
  18.   subroutines written in C or PASCAL. This is important for tedious
  19.   calculations like an optimization, where you have problems with the
  20.   accuracy. For this problems you can avoid the use of a super computer
  21.   which you feed with a C or PASCAL program for high precision numbers.
  22.  
  23. * All important functions of PASCAL and of a C library ("MATH.H") are
  24.   included. This includes the transcedental functions like sin(x), atan(x),
  25.   exp(x), log(x)....
  26.  
  27. * Programs written in C or C++, can use MFLOAT simply by replacing the 
  28.   double data format with the "MFLOAT" data format. The operators and the 
  29.   functions are overloaded for C++ and so only little changes of the source 
  30.   code are necessary. (You need BORLAND C 3.1).
  31.  
  32. * The subroutines are included in compiled form (object files) in a program.
  33.   Therfore you can use them for PASCAL, C and C++. The subroutines are
  34.   tested for BORLAND PASCAL 7.0, TURBO PASCAL 6.0 and BORLAND C 3.1.
  35.  
  36. * This package is limited to about 72 decimal digits. You may license the
  37.   source code. Then you can extend the precision for the algebraic 
  38.   calculations up to 10000 decimal digits (upper limit is due to the memory 
  39.   limit of one segment of the 8086) and for transcendental function up to 500 
  40.   digits.
  41.  
  42. * The internal data format is a binary floating point format. The length of 
  43.   the mantissa can be choosen by the user from one word to 15 words.
  44.   (1 word = 16 bits). This correspondeds to about 5, 10, 14, 19, 24,
  45.   29, 34, 39, 43, 48, 53, 58, 63, 67 and 72 digits of a decimal number.
  46.   1 word represents 4.81648 digits of a decimal number in average.
  47.   The exponent of two has the length of 16 bits. The largest representable
  48.   number is about 7.07E+9863, the smallest is about 7.07E-9865.
  49.  
  50. ***************************************************************************
  51.  
  52. The package includes following files:
  53.   READ.ME      - general information
  54.   MFLOAT.DOC   - describtion of the subroutines and mfloat numbers
  55.   MFLOATA.OBJ  - subroutines written in assembler
  56.   MFLOATB.OBJ  - subroutines written in C
  57.   MFLOAT.H     - header file for C
  58.   MFLOAT.CXX   - file for C++ : overloading functions and operators
  59.   PFLOAT.PAS   - unit: declarations of the external subroutines
  60.   PI.PAS       - example in PASCAL, calculation of pi
  61.   BESSEL.PAS   - example in PASCAL, calculation of the bessel function
  62.   POLAR.PAS    - example in PASCAL, conversion cartesian to polar coordinates
  63.   PI.C         - example in C, calculation of pi
  64.   BESSEL.C     - example in C, calculation of the bessel function
  65.   POLAR.C      - example in C, conversion cartesian to polar coordinates
  66.   PI.CPP       - example in C++, calculation of pi
  67.   BESSEL.CPP   - example in C++, calculation of the bessel function
  68.   POLAR.CPP    - example in C++, conversion cartesian to polar coordinates
  69.   PI.EXE       - example program executable
  70.  
  71. The source code option includes further files:
  72.   MFLOATA.ASM  - source code in assembler (basic subroutines)
  73.   MFLOATB.C    - source code in C (transcendental functions and extensions)
  74.   MFLOATC.ASM  - constants (include file of "MFLOATA.ASM")
  75.   CALCCON.PAS  - calculation of constants in "MFLOATC.ASM" in PASCAL source
  76.   CALCCON.EXE  - calculation of constants in "MFLOATC.ASM"
  77.  
  78. Registration:
  79.   see READ.ME
  80.  
  81. The authors of MFLOAT are interested in a feedback from the users of "MFLOAT".
  82. If you have problems or you find some bugs (the subroutines are checked very
  83. thorougly, but the probability, that a software of large dimension has no bug
  84. is not large in reality), please describe them and send us a mail per
  85. internet. You help us to remove this errors of the subroutines.
  86.  
  87. Many thanks for your help and much success with "MFLOAT".
  88.  
  89.  
  90. Kaufmann Friedrich & Mueller Walter
  91. Students at the Technical University of Graz
  92.  
  93. surface mail:
  94. Kaufmann Friedrich
  95. Schuetzenhofgasse 22
  96. A-8010 GRAZ
  97. AUSTRIA
  98. EUROPE
  99.  
  100. email address:
  101. fkauf@fstgds06.tu-graz.ac.at
  102.  
  103. ****************************************************************************
  104. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  105. ****************************************************************************
  106.  
  107.  
  108.         Multiprecision calculations with "MFLOAT"
  109.         =========================================
  110.  
  111. "MFLOAT" is a package of subroutines written in assembler for a 8086 micro-
  112. processor or compatibles. It is made for users having problems with the
  113. accuracy of standard IEEE arithmetic.
  114.  
  115. The precision of IEEE-754 standard floating point numbers is good enough
  116. for most technical and scientific applications. But sometimes a scientist
  117. or a mathematician has problems in his calculations because the algorithm or
  118. the mathematical problem is very ill-conditioned (for example finding the
  119. roots of a polynomial). As the floating point numbers are only appoximations
  120. of real or rational numbers, there is an additional error at every floating
  121. point operation.
  122.  
  123. There is no method known to represent every real number exactly in a
  124. binary format and therefore an approximation is necessary. A practical method
  125. is to extend the mantissa of the floating point number, which results in
  126. less error at every operation. But there are limitations because of higher
  127. memory needs for every floating point number and increased execution time.
  128. This is why a scientist having numerical problems of this kind is looking
  129. for very fast subroutines for multiprecision floating point numbers.
  130.  
  131. This package represents such a package of subroutines. It uses a binary
  132. data format because there are optimal instructions for the manipulation of
  133. binary numbers in the instruction set of the 8086. A binary data format has
  134. also the advantage, that the used memory is minimized. The only disadvantage
  135. is the higher effort that is necessary for the input and output of the
  136. numbers in ASCII-Format.
  137.  
  138. The main goal is a minimal computing time for very accurate numbers.
  139. Therefore no time is wasted for rounding (every result is truncated) and
  140. the handling of abnormal numbers as NANs or denormals. The only special
  141. number is zero. For saving time, there is no stack overflow check.
  142.  
  143. In case of an overflow or some other calculation error an error flag is set,
  144. which can be tested during or at the end of the calculation. Of course in
  145. such a case the numerical result is wrong. In case of an underflow the result
  146. is set to zero without any message.
  147.  
  148. A multiprecision number is consists of an array of words (1 word = 2 byte
  149. = 16 bit). The first word represents the exponent to the base of two. It is
  150. represented with an offset of 8000 Hex (= 32768) and ranges from 1 to
  151. FFFF Hex (= 65535). The special value 0 is used to represent the floating
  152. point number zero. In this case the rest of the array is arbitrary. The
  153. second word of the array represents the first word of mantissa. The mantissa
  154. is normalized for every floating point number and therefore the most
  155. significant bit of the first word of mantissa is always set. As there is no
  156. information in this bit, it is removed and at its place the sign bit is
  157. stored. The sign bit is one if the floating point number is negative, and
  158. zero if the floating point number is positive. The further words of
  159. the multiprecision number are words of mantissa. The length of mantissa
  160. ranging from one to a constant (15 in this release) can be choosen by the
  161. user.
  162.  
  163. Example : 8002 490F Hex  represents 3.141540528..
  164.  
  165. Every word of mantissa represents in average 4.8165 (= 16 * log10(2)) digits
  166. of a decimal number. Therefore you need 15 words of mantissa for a accuracy
  167. of 72 decimal digits and the whole number uses 16 words or 32 bytes.
  168.  
  169. If you have the source code, you can increase the accuracy. The upper limit
  170. of the accuracy is given by the used stack of the subroutines, which is
  171. about 3 times the memory of a "MFLOAT" number for the arithmetic operations
  172. addition, subtraction, multiplication and division. The stack is limited to
  173. one segment of the 8086 which is 65536 byte. A further limit for
  174. transcendental functions is 125 mantissa words (602 digits of a
  175. decimal number) for the series of the sine. If you extend this limit, you
  176. will get an error, for an argument near pi/4.
  177. "MFLOATA.OBJ" includes 19 "MFLOAT" constants. This constants are stored in
  178. the code segment, which is limited to 65536 bytes.
  179.  
  180. There is a very large range for the exponent. The greatest number, that can
  181. by represented is about +/-7.07E+9863, the smallest is about +/-7.07E-9865.
  182.  
  183. Now lets take a look inside the subroutines. The coprocessor 8087 or a
  184. compatible is not used because the IEEE format does not match the data
  185. format. There are subroutines for the algebraic operations addition,
  186. subtraction, multiplication and division.
  187.  
  188. Addition and subtraction are done by the same subroutine. The time for an
  189. addition (or subtraction) is proportional to the actual number of mantissa-
  190. words.
  191.  
  192. The implementation of the multiplication uses the instruction "MUL".
  193. which is much faster than shifted additions. The multiplication result
  194. is two times as large as its operands but only the half of the accuracy is
  195. needed. Therfore its a waste of time to calculate the second part of the
  196. mantissa. If exactly half the number of mantissawords are calculated and the
  197. rest is neglected, there is an accumulation of truncating errors, which can
  198. be very large. Therfore an additional extra word is added to the result during
  199. calculation, but for the final result this extra word is ignored. By this
  200. way half of the calculation time is saved and only small additional errors
  201. occure. The calculation time increases quadratic with the number of
  202. mantissawords.
  203.  
  204. The division is implemented using an algorithm similar to hand calculation.
  205. First an estimate of the first word of the result is made and the rest is
  206. calculated. If the rest is negativ, the estimation is corrected until it is
  207. ok. Using the rest the next word of mantissa is estimated and eventually
  208. corrected. This goes on until the whole mantissa of the result is calculated.
  209. For saving calculation time, the rest is truncated at every step by one word.
  210. Only at the first step an extra word is added (no truncation). Because of this
  211. method half of the calculation time is saved (same priciple as by the
  212. multiplication). The calculation time increases quadratic with the number of
  213. mantissawords.
  214.  
  215. The first transcendental function is the square root. There is a hand
  216. calculation method known similar to a division. The implemented algorithm
  217. uses this method. For saving calcultion time the same tricks are used as
  218. in the implementation of multiplication and division. The calculation time is
  219. smaller than those of a division. There are also iterative methods known to
  220. calculate the square root. These methods are much slower than the implemented
  221. method. (3 - 5 times slower).
  222.  
  223. Further trancendental functions are the trigonometric, hyperbolic functions,
  224. the exponential function and their invers functions. This functions can
  225. be derived from four series:
  226.  
  227. sin(x) = x - x^3/6 + x^5/120 - .....
  228.  
  229. sinh(x) = x + x^3/6 + x^5/120 + .....
  230.  
  231. arctan(x) = x - x^3/3 + x^5/5 - .....   |x| < 1
  232.  
  233. artanh(x) = x + x^3/3 + x^5/5 + .....   |x| < 1
  234.  
  235. These are the basic functions for further calculations. All other functions
  236. can be derived from them. Here are some examples:
  237.  
  238. cos(x) = sin(x + pi/2)
  239. tan(x) = sin(x) / cos(x)
  240. cot(x) = cos(x) / sin(x)
  241. cosh(x) = sqrt(1 + sqr(sinh(x))
  242. tanh(x) = sinh(x) / cosh(x)
  243. coth(x) = cosh(x) / sinh(x)
  244. exp(x) = sinh(x) + cosh(x)
  245. arcsin(x) = arctan(x / sqrt(1 - sqr(x)))
  246. arccos(x) = pi/2 - arcsin(x)
  247. arccot(x) = pi/2 - arctan(x)
  248. arsinh(x) = artanh(x / sqrt(1 + sqr(x)) = ln(x + sqrt(1 + sqr(x)))
  249. arcosh(x) = artanh(sqrt(sqr(x) - 1) / x) = ln(x + sqrt(sqr(x) - 1))
  250. arcoth(x) = artanh(1/x)
  251. ln(x) = 2 * artanh((x - 1) / (x + 1))
  252.  
  253. Some identities are used to save calculation time:
  254. sin(x + 2 * pi) = sin(x)
  255. exp(x + ln(2)) = 2 * exp(x)
  256. ....
  257. and so on.
  258.  
  259. Most of the time is used for the calculation of the series. Therefore it is
  260. important to write an optimal assembler subroutine for the whole series.
  261. This series could also be programmed directly by the former subroutines
  262. but using a specialized assembler subroutine the calculation time can be
  263. reduced to one third at even higher accuracy. The calculation time for a
  264. series increases with the third power of the number of mantissawords.
  265.  
  266.  
  267.  
  268. The subroutines are testet using the following compilers:
  269. Turbo Pascal 6.0, Borland Pascal 7.0, Borland C/C++ 3.1
  270.  
  271. You may use the library in combination with other compilers but in this
  272. case you may have to adabt the header files. Please note, that pascal
  273. calling conventions (upercase only, argument order, no underlines, ...)
  274. is implemented. All routines are 'far' with 'far' data pointers.
  275.  
  276.  
  277.  
  278. ===========================================================================
  279. ***************************************************************************
  280. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  281. ***************************************************************************
  282. ===========================================================================
  283.  
  284.       Using the subroutines for BORLAND C 3.1 :
  285.       =========================================
  286.  
  287. To make the functions available you have to include MFLOAT.H in your source
  288. and to specify MFLOATA.OBJ and MFLOATB.OBJ in your project file.
  289.  
  290. #define mfloatwords 16           /* maximum precision */
  291. typedef int mfloata[mfloatwords];/* mfloat array */
  292. #define mfloat_ int far *        /* pointer to array */
  293. #define mfloat mfloata           /* for use in 'C' programs */
  294.  
  295. (declared in "MFLOAT.H")
  296. This type declaration for multiprecision floating point numbers sets
  297. a maximal precision of 15 mantissawords or about 72 decimal digits. The
  298. maximum precision is an assembly time constant. In this release it is fixed
  299. to 15 mantissawords. It is possible to change the type declaration to reduce
  300. the number of mantissawords. This leads to a saving in memory needs and to
  301. higher execution speed at the cost of precision. It is also possible to adjust
  302. the precision dynamically using the function 'setmantissawords'.
  303.  
  304. The programmer must take care, that the stack is large enough, because
  305. there is no check of a stack overflow.
  306.  
  307. ****************************************************************************
  308. ++++++++++++++++++++     mfloat basic functions   ++++++++++++++++++++++++++
  309. ****************************************************************************
  310.  
  311. void  far pascal setmantissawords(int number);
  312.  
  313. This procedure sets the length of the mantissa of mfloat numbers. The number
  314. of word must be in the range of [1..15], otherwise is is bounded to these
  315. values. If the number of mantissawords is increased and the numbers of
  316. former calculations are used further, the additional words are not
  317. initialized to zero, their values are not known. The error of the number
  318. due to this further mantissa words is of course very small and is in the
  319. range of the accuracy of the former calculation. If it is important, that
  320. this further mantissa words are zero, the user has to initialize them (see
  321. data format). If the number of mantissa words is reduced, every number can
  322. be seen as truncated, although nothing is changed at this numbers. The number
  323. of mantissa words refers only to calculations.
  324.  
  325. ****************************************************************************
  326.  
  327. int   far pascal getmantissawords(void);
  328.  
  329. The function returns the actual number of mantissa words.
  330.  
  331. ****************************************************************************
  332.  
  333. void  far pascal reseterror(void);
  334.  
  335. The error flag is reset.
  336.  
  337. ****************************************************************************
  338.  
  339. char  far pascal geterror(void);
  340.  
  341. The error flag is returned. If there ocures a runtime error during the
  342. calculation, the error flag is set to true (=1).
  343.  
  344. ****************************************************************************
  345.  
  346. mfloat_  far pascal equm(mfloat_ a,                          /* a <-- b */
  347.              const mfloat_ b);
  348.  
  349. The value of b is copied to variable a with the defined accuracy. The value
  350. of b is not changed. The return value is a pointer to a.
  351.  
  352. ****************************************************************************
  353.  
  354. mfloat_  far pascal addm(mfloat_ a,                       /* a <-- a + b */
  355.              const mfloat_ b);
  356.  
  357. The value of b is added to the value a with the defined accuracy. The
  358. result is variable a. The value of b is not changed. The return value is a
  359. pointer to a.
  360.  
  361. ****************************************************************************
  362.  
  363. mfloat_  far pascal subm(mfloat_ a,                       /* a <-- a - b */
  364.              const mfloat_ b);
  365.  
  366. The value of b is subtracted from the value a with the defined accuracy.
  367. The result is variable a. The value of b is not changed. The return value is
  368. a pointer to a.
  369.  
  370. ****************************************************************************
  371.  
  372. mfloat_  far pascal multm(mfloat_ a,                      /* a <-- a * b */
  373.               const mfloat_ b);
  374.  
  375. The value of a is multiplied by value b with the defined accuracy. The result
  376. is variable a. The value of b is not changed. The return value is a pointer
  377. to a.
  378.  
  379. ****************************************************************************
  380.  
  381. mfloat_  far pascal divm(mfloat_ a,                       /* a <-- a / b */
  382.              const mfloat_ b);
  383.  
  384. The value of a is divided by the value b with the defined accuracy.
  385. The result is variable a. The value of b is not changed. The return value is
  386. a pointer to a.
  387.  
  388. ****************************************************************************
  389.  
  390. mfloat_  far pascal multi(mfloat_ a,                      /* a <-- a * b */
  391.                int b);
  392.  
  393. The value of a is multiplied by the integer value b with the defined
  394. accuracy. The result is variable a. The return value is a pointer to a.
  395.  
  396. ****************************************************************************
  397.  
  398. mfloat_  far pascal divi(mfloat_ a,                       /* a <-- a / b */
  399.              int b);
  400.  
  401. The value of a is divided by the integer value b with the defined accuracy.
  402. The result is variable a. The return value is a pointer to a.
  403.  
  404. ****************************************************************************
  405.  
  406. mfloat_  far pascal inversm(mfloat_ a);                   /* a <-- 1 / a */
  407.  
  408. The number one is divided by the number a with the defined accuracy.
  409. The result is variable a. The return value is a pointer to a.
  410.  
  411. ****************************************************************************
  412.  
  413. mfloat_  far pascal negm(mfloat_ a);                      /* a <--  - a  */
  414.  
  415. This procedure changes the sign of the number of variable a. The return value
  416. is a pointer to a.
  417.  
  418. ****************************************************************************
  419.  
  420. char  far pascal eqzero(const mfloat_ a);        /* eqzero <-- (a == 0)  */
  421.  
  422. This function checks, if the multiprecision number a is zero. The value of a
  423. is not changed. The return value is 0 or 1.
  424.  
  425. ****************************************************************************
  426.  
  427. char  far pascal gtzero(const mfloat_ a);         /* gtzero <-- (a > 0)  */
  428.  
  429. This function checks, if the multiprecision number a is greater than zero.
  430. The value of a is not changed. The return value is 0 or 1.
  431.  
  432. ****************************************************************************
  433.  
  434. char  far pascal gezero(const mfloat_ a);         /* gezero <-- (a >= 0) */
  435.  
  436. This function checks, if the multiprecision number a is greater than zero or
  437. equal zero. The value of a is not changed. The return value is 0 or 1.
  438.  
  439. ****************************************************************************
  440.  
  441. char  far pascal gtm(const mfloat_ a,                 /* gtm <-- (a > b) */
  442.              const mfloat_ b);
  443.  
  444. This function checks, if the multiprecision number a is greater than the
  445. multiprecision number b. The values of a and b are not changed. The return
  446. value is 0 or 1.
  447.  
  448. ****************************************************************************
  449.  
  450. char  far pascal eqm(const mfloat_ a,                /* eqm <-- (a == b) */
  451.              const mfloat_ b);
  452.  
  453. This function checks, if the multiprecision number a is equal to the
  454. multiprecision number b. The values of a and b are not changed. The return
  455. value is 0 or 1.
  456.  
  457. ****************************************************************************
  458.  
  459. mfloat_  far pascal getzerom(mfloat_ a);                     /* a <-- 0  */
  460.  
  461. This procedure copies the value zero = 0.0000.. to the variable a. The return
  462. value is a pointer to a.
  463.  
  464. ****************************************************************************
  465.  
  466. mfloat_  far pascal getonem(mfloat_ a);                       /* a <-- 1 */
  467.  
  468. This procedure copies the value one = 1.0000.. to the variable a. The return
  469. value is a pointer to a.
  470.  
  471. ****************************************************************************
  472.  
  473. mfloat_  far pascal getpim(mfloat_ a);                       /* a <-- pi */
  474.  
  475. This procedure copies the value pi = 3.14159.. to the variable a. The return
  476. value is a pointer to a.
  477.  
  478. ****************************************************************************
  479.  
  480. mfloat_  far pascal getln2m(mfloat_ a);                   /* a <-- ln(2) */
  481.  
  482. This procedure copies the value ln(2) = 0.69314.. to the variable a.
  483. The return value is a pointer to a.
  484.  
  485. ****************************************************************************
  486.  
  487. mfloat_  far pascal getln10m(mfloat_ a);                 /* a <-- ln(10) */
  488.  
  489. This procedure copies the value ln(10) = 2.3025.. to the variable a.
  490. The return value is a pointer to a.
  491.  
  492. ****************************************************************************
  493.  
  494. int  far pascal strtomf(mfloat_ a,                            /* a <-- b */
  495.             const char far * b);
  496.  
  497. This procedure converts the string b to the multiprecision number a. If the
  498. conversion is correct, the integer zero is returned. Otherwise the number
  499. of the character, where the first error occurs, is retured. For conversion
  500. the accuracy of the calculation is set to maximum and afterwards is set back
  501. to the original value. This function is not optimized for speed, but for
  502. accuracy, because it is not used very often normally.
  503.  
  504. ****************************************************************************
  505.  
  506. char far * far pascal mftoa(char far * str,                 /* str <-- a */
  507.                 const mfloat_ a,
  508.                 int len);
  509.  
  510. This procedure converts the multiprecision number a to a string str with
  511. the number of ASCII - characters, which is defined in len. If the number
  512. of characters is too less, the string is set to stars '*'. The format used
  513. for this conversion is '.32767F' (see mftostr). For conversion the
  514. accuracy of the calculation is set to maximum and afterwards is set back to
  515. the original value. This function is not optimized for speed, but for
  516. accuracy, because it is not used very often normally.
  517. The return value is a pointer to the character array a.
  518.  
  519. ****************************************************************************
  520.  
  521. char far * far pascal mftostr(char far * str,               /* str <-- a */
  522.                   const mfloat_ a,
  523.                   int len,
  524.                   const char far * format);
  525.  
  526. This procedure converts the multiprecision number a to a string str with
  527. a maximum number of ASCII - characters, which is defined in len. If the
  528. number of characters is too less, the string is set to stars '*'.
  529. The characters in the format - string have following meaning:
  530.  
  531. '.'   the point is shown ervery time, not significant zeros behind the point
  532.       are shown
  533. 'E'   scientific format with an exponent of ten
  534. 'F'   fixpoint format, if the number is too large, a scientific format
  535.       with an exponent of ten is used.
  536. 'G'   fixpoint format, if the number is too large, or is less than 0.0001,
  537.       a scientific format is used.
  538. '1..' an integer, which declares the number of digits behind the point
  539.       for fixpoint and scientific format.
  540.  
  541. The order of the characters (but not the digits for the precision) in the
  542. format string is arbitrary. It is also allowed to use lower case characters
  543. for the format, but in this case the letter for the exponent of ten is
  544. also the lower case letter 'e'. For conversion the accuracy of the
  545. calculation is set to maximum and afterwards is set back to the original
  546. value. This function is not optimized for speed, but for accuracy, because
  547. it is not used very often normally.
  548. The return value is a pointer to the character array a.
  549.  
  550. ****************************************************************************
  551.  
  552. double far pascal mftod(const mfloat_ b);                 /* mftod <-- b */
  553.  
  554. This function converts a multiprecision number to a double number.
  555. The value of b is not changed.
  556.  
  557. ****************************************************************************
  558.  
  559. long double far pascal mftold(const mfloat_ b);          /* mftold <-- b */
  560.  
  561. This function converts a multiprecision number to a long double number.
  562. The value of b is not changed.
  563.  
  564. ****************************************************************************
  565.  
  566. mfloat_  far pascal dtomf(mfloat_ a,                          /* a <-- b */
  567.               double b);
  568.  
  569. This procedure converts a double precision number to a multiprecision number.
  570. The return value is a pointer to a.
  571.  
  572. ****************************************************************************
  573.  
  574. mfloat_  far pascal ldtomf(mfloat_ a,                         /* a <-- b */
  575.                long double b);
  576.  
  577. This procedure converts a long double number to a multiprecision number.
  578. The return value is a pointer to a.
  579.  
  580. ****************************************************************************
  581. +++++++++++++++     standart functions (Borland C: MATH.H)   +++++++++++++++
  582. ****************************************************************************
  583.  
  584. mfloat_  far pascal acosm(mfloat_ a);                   /* a <-- acos(a) */
  585.  
  586. This procedure calculates the invers cos(x) function of the argument a. The
  587. result is the variable a.
  588. -1 <= a <= 1        range of the argument, otherwise the error flag is set
  589.  0 <= a <= pi       range of the result
  590. The return value is a pointer to a.
  591.  
  592. ****************************************************************************
  593.  
  594. mfloat_  far pascal asinm(mfloat_ a);                   /* a <-- asin(a) */
  595.  
  596. This procedure calculates the invers sin(x) function of the argument a. The
  597. result is the variable a.
  598. -1 <= a <= 1        range of the argument, otherwise the error flag is set
  599. -pi/2 <= a <= pi/2  range of the result
  600. The return value is a pointer to a.
  601.  
  602. ****************************************************************************
  603.  
  604. mfloat_  far pascal atanm(mfloat_ a);                   /* a <-- atan(a) */
  605.  
  606. This procedure calculates the invers tan(x) function of the argument a. The
  607. result is the variable a.
  608. -pi/2 < a < pi/2       range of the result
  609. The return value is a pointer to a.
  610.  
  611. ****************************************************************************
  612.  
  613. mfloat_  far pascal atan2m(mfloat_ a,                /* a <-- atan2(a,b) */
  614.                const mfloat_ b);
  615.  
  616. This function calulates the extended arctan function phi = arctan(a/b). If
  617. a is zero or less than zero, the correct value phi is retured, which is
  618. in the intervall (-pi, pi]:
  619. a = r * sin(phi)
  620. b = r * cos(phi)
  621. r = sqrt(sqr(a) + sqr(b)) >= 0
  622.  
  623. If a = 0 and b = 0, the value phi is not unique. The error flag is set and
  624. the number zero is returned. The return variable is the variable a. Variable
  625. b is not changed.
  626. The return value is a pointer to a.
  627.  
  628. ****************************************************************************
  629.  
  630. mfloat_  far pascal ceilm(mfloat_ a);                   /* a <-- ceil(a) */
  631.  
  632. The least integer number, which is greater than or equal the mfloat number
  633. a, is calculated and returned to the variable a.
  634. The return value is a pointer to a.
  635.  
  636. ****************************************************************************
  637.  
  638. mfloat_  far pascal cosm(mfloat_ a);                     /* a <-- cos(a) */
  639.  
  640. This procedure calculates the cos(x) function of the argument a. The result
  641. is the variable a.
  642. The return value is a pointer to a.
  643.  
  644. ****************************************************************************
  645.  
  646. mfloat_  far pascal coshm(mfloat_ a);                   /* a <-- cosh(a) */
  647.  
  648. This procedure calculates the cosh(x) function of the argument a. The result
  649. is the variable a.
  650. -22712 < a < 22712    range of the argument without an overflow
  651.      a >= 1       range of the result
  652. The return value is a pointer to a.
  653.  
  654. ****************************************************************************
  655.  
  656. mfloat_  far pascal expm(mfloat_ a);                     /* a <-- exp(a) */
  657.  
  658. This procedure calculates exponential function exp(x) of the argument a. The
  659. result is the variable a.
  660.  
  661.  a < 22712   range of the argument without an overflow
  662.  a < -22712  there is an underflow result = 0
  663.  a >= 0      range of the result
  664. The return value is a pointer to a.
  665.  
  666. ****************************************************************************
  667.  
  668. mfloat_  far pascal fabsm(mfloat_ a);                   /* a <-- fabs(a) */
  669.  
  670. This procedure returns the absolute value of the number a. The result is
  671. the variable a.
  672. The return value is a pointer to a.
  673.  
  674. ****************************************************************************
  675.  
  676. mfloat_  far pascal floorm(mfloat_ a);                 /* a <-- floor(a) */
  677.  
  678. The greatest integer number, which is less than or equal the mfloat number
  679. a, is calculated and returned to the variable a.
  680. The return value is a pointer to a.
  681.  
  682. ****************************************************************************
  683.  
  684. mfloat_  far pascal fmodm(mfloat_ a,                  /* a <-- fmod(a,b) */
  685.               const mfloat_ b);
  686.  
  687. This is a floating point modul function. The number a is divided by b, so
  688. that the result is an integer. The rest of this division is returned by the
  689. variable a. The variable b is not changed.
  690.  
  691. a = k * b + c;   k is an integer;   0 <= c / b < 1;  c is assigned to a
  692.  
  693. If the variable b is zero, the error flag is set. For negative b, the result
  694. is negativ or zero. For very large k, which can not be represented as an
  695. mfloat number exactly, the result is zero.
  696. The return value is a pointer to a.
  697.  
  698. ****************************************************************************
  699.  
  700. mfloat_  far pascal frexpm(mfloat_ a,                /* a <-- frexp(a,b) */
  701.                int * b);
  702.  
  703. The mfloat number a is splitted in the exponent and in the mantissa:
  704.  
  705. a = c * (2**b);  0.5 <= |c| < 1 if a <> 0;  c is assigned to a
  706.  
  707. This procedure returns two results. IF the mfloat number is zero, it remains
  708. zero and b is set to zero.
  709. The return value is a pointer to a.
  710.  
  711. ****************************************************************************
  712.  
  713. mfloat_  far pascal hypotm(mfloat_ a,                /* a <-- hypot(a,b) */
  714.                const mfloat_ b);
  715.  
  716. This function is used for the calculation of the absolute value of a complex
  717. number:
  718.  
  719. c = sqrt(sqr(a) + sqr(b)); c is assigned to a; b remains unchanged.
  720. |a + i*b| = hypot(a,b);
  721. This procedure calculates the correct result also if sqr(a) can not be
  722. represented as a mfloat number, because it is too great or too small.
  723. The return value is a pointer to a.
  724.  
  725. ****************************************************************************
  726.  
  727. mfloat_  far pascal ldexpm(mfloat_ a,                /* a <-- ldexp(a,b) */
  728.                int b);
  729.  
  730. The mfloat number a is multiplied by the power of two 2**b.
  731.  
  732. a <-- a * (2**b);
  733.  
  734. The calculation is very fast because only the exponent of two is manipulated.
  735. The return value is a pointer to a.
  736.  
  737. ****************************************************************************
  738.  
  739. mfloat_  far pascal logm(mfloat_ a);                     /* a <-- log(a) */
  740.  
  741. This procedure calculates the natural logarithm ln(x) of the argument a. The
  742. result is the variable a.
  743.  
  744. a > 0                  range of the argument
  745. -22712 <= a <= 22712   range of the result
  746. The return value is a pointer to a.
  747.  
  748. ****************************************************************************
  749.  
  750. mfloat_  far pascal log10m(mfloat_ a);                 /* a <-- log10(a) */
  751.  
  752. This procedure calculates the decade logarithm log(x) of the argument a. The
  753. result is the variable a.
  754.  
  755. a > 0                  range of the argument
  756. -9864 <= a <= 9863     range of the result
  757. The return value is a pointer to a.
  758.  
  759. ****************************************************************************
  760.  
  761. mfloat_  far pascal modfm(mfloat_ a,                  /* a,b <-- modf(a) */
  762.               mfloat_ b);
  763.  
  764. This procedure splits the number a in an integer part and in an fractional
  765. part:
  766.  
  767. a = k + b;    k is an integer;    0 <= b < 1;    k is assigned to a
  768. The return value is a pointer to a.
  769.  
  770. ****************************************************************************
  771.  
  772. mfloat_  far pascal powm(mfloat_ a,                    /* a <-- pow(a,b) */
  773.              const mfloat_ b);
  774.  
  775. This procedure calculates the power of the basis a by the exponent b.
  776. The result is assigned to the variable a, the variable b is not changed.
  777. If the exponent is an integer, the calculation is made by repeated
  778. multiplications. The number of multiplications is proportional to
  779. log(|b|). If b has a fractional, the calculation is made by logarithm. In
  780. this case, the base a must be positiv, otherwise the error flag is set,
  781. because the result is not real. If the number b is very large so that
  782. no fractional can be represented by the mfloat format and the base
  783. is negativ, the error flag is set too. If a = 0 and b <= 0 the result is
  784. not defined and the error flag is set.
  785. The return value is a pointer to a.
  786.  
  787. ****************************************************************************
  788.  
  789. mfloat_  far pascal pow10m(mfloat_ a,                     /* a <-- 10**b */
  790.                int b);
  791.  
  792. The power of the basis ten is build with the integer exponent b. The result
  793. is assigned to a.
  794.  
  795. a <-- 10**b
  796.  
  797. if b < -9864 then a = 0
  798. if b > 9863  then there is an overflow -> error flag is set
  799. The return value is a pointer to a.
  800.  
  801. ****************************************************************************
  802.  
  803. mfloat_  far pascal sinm(mfloat_ a);                     /* a <-- sin(a) */
  804.  
  805. This procedure calculates the sin(x) function of the argument a. The result
  806. is the variable a.
  807. The return value is a pointer to a.
  808.  
  809. ****************************************************************************
  810.  
  811. mfloat_  far pascal sinhm(mfloat_ a);                   /* a <-- sinh(a) */
  812.  
  813. This procedure calculates the sinh(x) function of the argument a. The result
  814. is the variable a.
  815. -22712 < a < 22712   range of the argument without an overflow
  816. The return value is a pointer to a.
  817.  
  818. ****************************************************************************
  819.  
  820. mfloat_  far pascal sqrtm(mfloat_ a);                   /* a <-- sqrt(a) */
  821.  
  822. This procedure calculates the square root of the number a. The result is
  823. variable a. If a is negativ, the error flag is set.
  824. The return value is a pointer to a.
  825.  
  826. ****************************************************************************
  827.  
  828. mfloat_  far pascal tanm(mfloat_ a);                     /* a <-- tan(a) */
  829.  
  830. This procedure calculates the tan(x) function of the argument a. The result
  831. is the variable a. It is faster than the calculation of sin(a) / cos(a).
  832. The return value is a pointer to a.
  833.  
  834. ****************************************************************************
  835.  
  836. mfloat_  far pascal tanhm(mfloat_ a);                   /* a <-- tanh(a) */
  837.  
  838. This procedure calculates the tanh(x) function of the argument a. The result
  839. is the variable a.
  840. The return value is a pointer to a.
  841.  
  842. ****************************************************************************
  843. ++++++++++++++++++    extended standart functions   ++++++++++++++++++++++++
  844. ****************************************************************************
  845.  
  846. mfloat_  far pascal acoshm(mfloat_ a);                 /* a <-- acosh(a) */
  847.  
  848. This procedure calculates the invers cosh(x) function of the argument a. The
  849. result is the variable a.
  850. 1 <= a          range of the argument
  851. 0 <= a < 22712  range of the result
  852. The return value is a pointer to a.
  853.  
  854. ****************************************************************************
  855.  
  856. mfloat_  far pascal acotm(mfloat_ a);                   /* a <-- acot(a) */
  857.  
  858. This procedure calculates the invers cot(x) function of the argument a. The
  859. result is the variable a.
  860. 0 < a < pi       range of the result
  861. The return value is a pointer to a.
  862.  
  863. ****************************************************************************
  864.  
  865. mfloat_  far pascal acothm(mfloat_ a);                 /* a <-- acoth(a) */
  866.  
  867. This procedure calculates the invers coth(x) function of the argument a. The
  868. result is the variable a.
  869.  -1 > a   or   1 < a    range of the argument
  870. The return value is a pointer to a.
  871.  
  872. ****************************************************************************
  873.  
  874. mfloat_  far pascal asinhm(mfloat_ a);                 /* a <-- asinh(a) */
  875.  
  876. This procedure calculates the invers sinh(x) function of the argument a. The
  877. result is the variable a.
  878. The return value is a pointer to a.
  879.  
  880. ****************************************************************************
  881.  
  882. mfloat_  far pascal atanhm(mfloat_ a);                 /* a <-- atanh(a) */
  883.  
  884. This procedure calculates the invers tanh(x) function of the argument a. The
  885. result is the variable a.
  886. -1 < a < 1  range of the argument
  887. The return value is a pointer to a.
  888.  
  889. ****************************************************************************
  890.  
  891. mfloat_  far pascal cossinm(mfloat_ a,                  /* a <-- cos(a), */
  892.                 mfloat_ b);                 /* b <-- sin(a)  */
  893.  
  894. This procedure calculates the sin(x) and cos(x) function by only one
  895. evulation of the series for the sine function. This procedure is therefore
  896. faster than two individual calculations. It is recommended for example
  897. for the calculation of the exponential function of complex arguments.
  898. The return value is a pointer to a.
  899.  
  900. ****************************************************************************
  901.  
  902. mfloat_  far pascal cotm(mfloat_ a);                     /* a <-- cot(a) */
  903.  
  904. This procedure calculates the cot(x) function of the argument a. The result
  905. is the variable a. It is faster than the calculation of cos(a) / sin(a).
  906. The return value is a pointer to a.
  907.  
  908. ****************************************************************************
  909.  
  910. mfloat_  far pascal cothm(mfloat_ a);                   /* a <-- coth(a) */
  911.  
  912. This procedure calculates the coth(x) function of the argument a. The result
  913. is the variable a.
  914. range of the result: a >=1 or a <= -1
  915. The return value is a pointer to a.
  916.  
  917. ****************************************************************************
  918.  
  919. mfloat_  far pascal exp10m(mfloat_ a);                  /* a <-- 10 ** a */
  920.  
  921. This procedure calculates the power of 10 of the argument a. The
  922. result is the variable a.
  923.  a < 9863    range of the argument without an overflow
  924.  a < -9864   there is an underflow, result = 0
  925.  a >= 0      range of the result
  926. The return value is a pointer to a.
  927.  
  928. ****************************************************************************
  929.  
  930. mfloat_  far pascal sqrm(mfloat_ a);                     /* a <-- sqr(a) */
  931.  
  932. This procedure calculates the sqare of a. The result is variabel a.
  933. sqr(a) = a * a
  934. The return value is a pointer to a.
  935.  
  936. ****************************************************************************
  937.  
  938. mfloat_  far pascal truncm(mfloat_ a);                 /* a <-- trunc(a) */
  939.  
  940. This procedure returns the integer part of the number a (truncation). The
  941. result is variable a. The sign has no influence and remains unchanged.
  942. The return value is a pointer to a.
  943.  
  944.  
  945.  
  946.  
  947.  
  948. ===========================================================================
  949. ***************************************************************************
  950. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  951. ***************************************************************************
  952. ===========================================================================
  953.  
  954.   Using the subroutines for Borland C++ 3.1 :
  955.   ===========================================
  956.  
  957. Using the mfloat library with C++ all of the mathematical operators are
  958. overloaded to allow 'normal' calculations with the mfloat class.
  959.  
  960. To make the functions available you have to include MFLOAT.CXX (includes
  961. MFLOAT.H) in your source and to specify MFLOATA.OBJ and MFLOATB.OBJ
  962. in your project file.
  963.  
  964. All operators available for the type 'double' are also defined for 'mfloat',
  965. So it is very easy to convert a program using double variables to use the
  966. mfloat library. Please note that it is not allowed to assign 'mfloat'
  967. variables to other numeric types. To solve this problem please use 'mftod'
  968. and 'mftold' (see the example programs!).
  969.  
  970. overloaded operators:
  971.   =
  972.   ++ (prefix, suffix)     -- (prefix, suffix)
  973.   + (unary, binary)       -  (unary, binary)
  974.   *     /
  975.   +=    -=    *=    /=
  976.   <     >     <=    >=    ==    !=
  977.   <<  (output to stream)  >>  (input from stream)
  978.  
  979. The following functions are overloaded mfloat equivalents of the
  980. functions defined in 'math.h':
  981.  
  982. acos, asin, atan, atan2, ceil, cos, cosh, exp, fabs, floor, fmod, frexp,
  983. hypot, ldexp, log, log10, modf, pow, sin, sinh, sqrt, tan, tanh.
  984.  
  985. Instead of 'atof' use 'atofm' (overloading problem).
  986. Instead of 'pow10' use 'pow10m' (overloading problem).
  987.  
  988. The following functions are extensions to math.h:
  989.  
  990. acosh, asinh, atanh, acoth, cot, coth, exp10, sqr, trunc.
  991.  
  992. Further functions are available:
  993.  
  994. void   pascal setmantissawords(int number);
  995. int    pascal getmantissawords(void);
  996. void   pascal reseterror(void);
  997. char   pascal geterror(void);
  998. int    pascal strtomf(mfloat &a, const char *s);
  999. char * pascal mftoa(char *str, const mfloat &a, int len);
  1000. char * pascal mftostr(char *str, const mfloat &a, int len,
  1001.               const char * format);
  1002. double pascal mftod(const mfloat &b);
  1003. long double pascal mftold(const mfloat &b);
  1004.  
  1005. All functions are mapped to the functions defined in MFLOAT.H. So please
  1006. take a look in description for using the functions with Borland C above.
  1007.  
  1008.  
  1009.  
  1010.  
  1011. ===========================================================================
  1012. ***************************************************************************
  1013. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  1014. ***************************************************************************
  1015. ===========================================================================
  1016.  
  1017.   Using the subroutines for TURBO PASCAL 6.0, 7.0 or BORLAND PASCAL 7.0 :
  1018.   =======================================================================
  1019.  
  1020. CONST MfloatWords = 16;
  1021. TYPE mfloat = ARRAY[0..MfloatWords-1] OF integer;
  1022.  
  1023. (declared in "PFLOAT.PAS")
  1024.  
  1025. This type declaration for multiprecision floating point numbers sets
  1026. a maximal precision of 15 mantissawords or about 72 decimal digits. The
  1027. maximum precision is an assembly time constant. In this release it is fixed
  1028. to 15 mantissawords. It is possible to change the type declaration to reduce
  1029. the number of mantissawords. This leads to a saving in memory needs and to
  1030. higher execution speed at the cost of precision. It is also possible to adjust
  1031. the precision dynamically using the function 'setmantissawords'.
  1032.  
  1033. The programmer must take care, that the stack is large enough, because
  1034. there is no check of a stack overflow.
  1035.  
  1036. ****************************************************************************
  1037. ++++++++++++++++++++     mfloat basic functions   ++++++++++++++++++++++++++
  1038. ****************************************************************************
  1039.  
  1040. void  far pascal setmantissawords(int number);
  1041.  
  1042. ****************************************************************************
  1043. ++++++++++++++++++++     mfloat basic functions   ++++++++++++++++++++++++++
  1044. ****************************************************************************
  1045.  
  1046. PROCEDURE Setmantissawords(number : integer);
  1047.  
  1048. This procedure sets the length of the mantissa of mfloat numbers. The number
  1049. of word must be in the range of [1..15], otherwise is is bounded to these
  1050. values. If the number of mantissawords is increased and the numbers of
  1051. former calculations are used further, the additional words are not
  1052. initialized to zero, their values are not known. The error of the number
  1053. due to this further mantissa words is of course very small and is in the
  1054. range of the accuracy of the former calculation. If it is important, that
  1055. this further mantissa words are zero, the user has to initialize them (see
  1056. data format). If the number of mantissa words is reduced, every number can
  1057. be seen as truncated, although nothing is changed at this numbers. The number
  1058. of mantissa words refers only to calculations.
  1059.  
  1060. ****************************************************************************
  1061.  
  1062. FUNCTION GetMantissawords : integer;
  1063.  
  1064. The function returns the actual number of mantissa words.
  1065.  
  1066. ****************************************************************************
  1067.  
  1068. PROCEDURE ResetError;
  1069.  
  1070. The error flag is reset.
  1071.  
  1072. ****************************************************************************
  1073.  
  1074. FUNCTION GetError : boolean;
  1075.  
  1076. The error flag is returned. If there ocures a runtime error during the
  1077. calculation, the error flag is set to true.
  1078.  
  1079. ****************************************************************************
  1080.  
  1081. PROCEDURE equm(VAR a, b : mfloat);                   { *** a <-- b *** }
  1082.  
  1083. The value of b is copied to variable a with the defined accuracy. The value
  1084. of b is not changed.
  1085.  
  1086. ****************************************************************************
  1087.  
  1088. PROCEDURE addm(VAR a, b : mfloat);                { *** a <-- a + b *** }
  1089.  
  1090. The value of b is added to the value a with the defined accuracy. The
  1091. result is variable a. The value of b is not changed.
  1092.  
  1093. ****************************************************************************
  1094.  
  1095. PROCEDURE subm(VAR a, b : mfloat);                { *** a <-- a - b *** }
  1096.  
  1097. The value of b is subtracted from the value a with the defined accuracy.
  1098. The result is variable a. The value of b is not changed.
  1099.  
  1100. ****************************************************************************
  1101.  
  1102. PROCEDURE multm(VAR a, b : mfloat);               { *** a <-- a * b *** }
  1103.  
  1104. The value of a is multiplied by value b with the defined accuracy.
  1105. The result is variable a. The value of b is not changed.
  1106.  
  1107. ****************************************************************************
  1108.  
  1109. PROCEDURE divm(VAR a, b : mfloat);                { *** a <-- a / b *** }
  1110.  
  1111. The value of a is divided by the value b with the defined accuracy.
  1112. The result is variable a. The value of b is not changed.
  1113.  
  1114. ****************************************************************************
  1115.  
  1116. PROCEDURE multi(VAR a : mfloat; b : integer);     { *** a <-- a * b *** }
  1117.  
  1118. The value of a is multiplied by the integer value b with the defined
  1119. accuracy. The result is variable a.
  1120.  
  1121. ****************************************************************************
  1122.  
  1123. PROCEDURE divi(VAR a : mfloat; b : integer);      { *** a <-- a / b *** }
  1124.  
  1125. The value of a is divided by the integer value b with the defined accuracy.
  1126. The result is variable a.
  1127.  
  1128. ****************************************************************************
  1129.  
  1130. PROCEDURE inversm(VAR a : mfloat);              { *** a <-- 1 / a *** }
  1131.  
  1132. The number one is divided by the number a with the defined accuracy.
  1133. The result is variable a.
  1134.  
  1135. ****************************************************************************
  1136.  
  1137. PROCEDURE negm(VAR a : mfloat);                 { *** a <--  - a *** }
  1138.  
  1139. This procedure changes the sign of the number of variable a.
  1140.  
  1141. ****************************************************************************
  1142.  
  1143. FUNCTION eqZero(VAR a : mfloat) : boolean;    { *** eqZero <-- a = 0  *** }
  1144.  
  1145. This function checks, if the multiprecision number a zero. The value of a
  1146. is not changed.
  1147.  
  1148. ****************************************************************************
  1149.  
  1150. FUNCTION gtZero(VAR a : mfloat) : boolean;    { *** gtZero <-- a > 0 *** }
  1151.  
  1152. This function checks, if the multiprecision number a is greater than zero.
  1153. The value of a is not changed.
  1154.  
  1155. ****************************************************************************
  1156.  
  1157. FUNCTION geZero(VAR a : mfloat) : boolean;    { *** geZero <-- a >= 0 *** }
  1158.  
  1159. This function checks, if the multiprecision number a is greater than zero or
  1160. equal zero. The value of a is not changed.
  1161.  
  1162. ****************************************************************************
  1163.  
  1164. FUNCTION gtm(VAR a, b : mfloat) : boolean;     { *** gtm <-- a > b *** }
  1165.  
  1166. This function checks, if the multiprecision number a is greater than the
  1167. multiprecision number b. The values of a and b are not changed.
  1168.  
  1169. ****************************************************************************
  1170.  
  1171. FUNCTION eqm(VAR a, b : mfloat) : boolean;     { *** eqm <-- a = b *** }
  1172.  
  1173. This function checks, if the multiprecision number a is equal to the
  1174. multiprecision number b. The values of a and b are not changed.
  1175.  
  1176. ****************************************************************************
  1177.  
  1178. PROCEDURE GetZerom(VAR a : mfloat);             { *** a <- 0 *** }
  1179.  
  1180. This procedure copies the value zero = 0.0000.. to the variable a.
  1181.  
  1182. ****************************************************************************
  1183.  
  1184. PROCEDURE GetOnem(VAR a : mfloat);              { *** a <- 1 *** }
  1185.  
  1186. This procedure copies the value one = 1.0000.. to the variable a.
  1187.  
  1188. ****************************************************************************
  1189.  
  1190. PROCEDURE GetPim(VAR a : mfloat);               { *** a <- pi *** }
  1191.  
  1192. This procedure copies the value pi = 3.14159.. to the variable a.
  1193.  
  1194. ****************************************************************************
  1195.  
  1196. PROCEDURE GetLn2m(VAR a : mfloat);              { *** a <- ln(2) *** }
  1197.  
  1198. This procedure copies the value ln(2) = 0.69314.. to the variable a.
  1199.  
  1200. ****************************************************************************
  1201.  
  1202. PROCEDURE GetLn10m(VAR a : mfloat);             { *** a <- ln(10) *** }
  1203.  
  1204. This procedure copies the value ln(10) = 2.3025.. to the variable a.
  1205.  
  1206. ****************************************************************************
  1207.  
  1208. FUNCTION strtomf(VAR a : mfloat;               { *** a <-- string *** }
  1209.              b : string) : integer;
  1210.  
  1211. This procedure converts the string b to the multiprecision number a. If the
  1212. conversion is correct, the integer zero is returned. Otherwise the number
  1213. of the character, where the first error occurs, is retured. For conversion
  1214. the accuracy of the calculation is set to maximum and afterwards is set back
  1215. to the original value. This function is not optimized for speed, but for
  1216. accuracy, because it is not used very often normally.
  1217.  
  1218. ****************************************************************************
  1219.  
  1220. FUNCTION  mftoa(VAR a   : mfloat;               { *** string <-- a *** }
  1221.             len : integer) : string;
  1222.  
  1223. This procedure converts the multiprecision number a to a string str with
  1224. the number of ASCII - characters, which is defined in len. If the number
  1225. of characters is too less, the string is set to stars '*'. The format used
  1226. for this conversion is '.32767F' (see FUNCTION mftostr). For conversion the
  1227. accuracy of the calculation is set to maximum and afterwards is set back to
  1228. the original value. This function is not optimized for speed, but for
  1229. accuracy, because it is not used very often normally.
  1230.  
  1231. ****************************************************************************
  1232.  
  1233. FUNCTION  mftostr(VAR a      : mfloat;           { *** string <-- a *** }
  1234.               len    : integer;
  1235.               format : string) : string;
  1236.  
  1237. This procedure converts the multiprecision number a to a string str with
  1238. a maximum number of ASCII - characters, which is defined in len. If the
  1239. number of characters is too less, the string is set to stars '*'.
  1240. The characters in the format - string have following meaning:
  1241.  
  1242. '.'   the point is shown ervery time, not significant zeros behind the point
  1243.       are shown
  1244. 'E'   scientific format with an exponent of ten
  1245. 'F'   fixpoint format, if the number is too large, a scientific format
  1246.       with an exponent of ten is used
  1247. 'G'   fixpoint format, if the number is too large, or less than 0.0001,
  1248.       a scientific format is used.
  1249. '1..' an integer, which declares the number of digits behind the point
  1250.       for fixpoint and scientific format
  1251.  
  1252. The order of the characters (but not the digits for the precision) in the
  1253. format string is arbitrary. It is also allowed to use lower case characters
  1254. for the format, but in this case the letter for the exponent of ten is
  1255. also the lower case letter 'e'. For conversion the accuracy of the
  1256. calculation is set to maximum and afterwards is set back to the original
  1257. value. This function is not optimized for speed, but for accuracy, because
  1258. it is not used very often normally.
  1259.  
  1260. ****************************************************************************
  1261.  
  1262. FUNCTION MfToD(VAR a : mfloat) : double;           { *** MfToD <- a *** }
  1263.  
  1264. This function converts a multiprecision number to a double precision number.
  1265. The value of a is not changed.
  1266.  
  1267. ****************************************************************************
  1268.  
  1269. FUNCTION MfToLd(VAR a : mfloat) : extended;       { *** MfToLd <- a *** }
  1270.  
  1271. This function converts a multiprecision number to an extended number. The
  1272. value of a is not changed.
  1273.  
  1274. ****************************************************************************
  1275.  
  1276. PROCEDURE DToMf(VAR a : mfloat; b : double);  { *** a <- b *** }
  1277.  
  1278. This procedure converts a double precision number to a multiprecision number.
  1279.  
  1280. ****************************************************************************
  1281.  
  1282. PROCEDURE LdToMf(VAR a : mfloat; b : extended);{ *** a <- b *** }
  1283.  
  1284. This procedure converts an extended number to a multiprecision number.
  1285.  
  1286. ****************************************************************************
  1287. ++++++++++++++++     standart functions (Borland C: MATH.H)   ++++++++++++++
  1288. ****************************************************************************
  1289.  
  1290. PROCEDURE acosm(VAR a : mfloat);                 { *** a <- arccos(a) *** }
  1291.  
  1292. This procedure calculates the invers cos(x) function of the argument a. The
  1293. result is the variable a.
  1294. -1 <= a <= 1        range of the argument, otherwise the error flag is set
  1295.  0 <= a <= pi       range of the result
  1296.  
  1297. ****************************************************************************
  1298.  
  1299. PROCEDURE asinm(VAR a : mfloat);                 { *** a <- arcsin(a) *** }
  1300.  
  1301. This procedure calculates the invers sin(x) function of the argument a. The
  1302. result is the variable a.
  1303. -1 <= a <= 1        range of the argument, otherwise the error flag is set
  1304. -pi/2 <= a <= pi/2  range of the result
  1305.  
  1306. ****************************************************************************
  1307.  
  1308. PROCEDURE atanm(VAR a : mfloat);                 { *** a <- arctan(a) *** }
  1309.  
  1310. This procedure calculates the invers tan(x) function of the argument a. The
  1311. result is the variable a.
  1312. -pi/2 < a < pi/2       range of the result
  1313.  
  1314. ****************************************************************************
  1315.  
  1316. PROCEDURE atan2m(VAR a, b : mfloat);           { *** a <- arg(b + ia) *** }
  1317.  
  1318. This function calulates the extended arctan function phi = arctan(a/b). If
  1319. a is zero or less than zero, the correct value phi is retured, which is
  1320. in the intervall (-pi, pi]:
  1321. a = r * sin(phi)
  1322. b = r * cos(phi)
  1323. r = sqrt(sqr(a) + sqr(b)) >= 0
  1324.  
  1325. If a = 0 and b = 0, the value phi is not unique. The error flag is set and
  1326. the number zero is returned. The return variable is the variable a. Variable
  1327. b is not changed.
  1328.  
  1329. ****************************************************************************
  1330.  
  1331. PROCEDURE ceilm(VAR a : mfloat);                 { *** a <-- ceil(a) *** }
  1332.  
  1333. The least integer number, which is greater than or equal the mfloat number
  1334. a, is calculated and returned to the variable a.
  1335.  
  1336. ****************************************************************************
  1337.  
  1338. PROCEDURE cosm(VAR a : mfloat);                    { *** a <- cos(a) *** }
  1339.  
  1340. This procedure calculates the cos(x) function of the argument a. The result
  1341. is the variable a.
  1342.  
  1343. ****************************************************************************
  1344.  
  1345. PROCEDURE coshm(VAR a : mfloat);                  { *** a <- cosh(a) *** }
  1346.  
  1347. This procedure calculates the cosh(x) function of the argument a. The result
  1348. is the variable a.
  1349. -22712 < a < 22712    range of the argument without an overflow
  1350.      a >= 1       range of the result
  1351.  
  1352. ****************************************************************************
  1353.  
  1354. PROCEDURE expm(VAR a : mfloat);                    { *** a <- exp(a) *** }
  1355.  
  1356. This procedure calculates exponential function exp(x) of the argument a. The
  1357. result is the variable a.
  1358.  
  1359.  a < 22712   range of the argument without an overflow
  1360.  a < -22712  there is an underflow result = 0
  1361.  a >= 0      range of the result
  1362.  
  1363. ****************************************************************************
  1364.  
  1365. PROCEDURE fabsm(VAR a : mfloat);                 { *** a <-- fabs(a) *** }
  1366.  
  1367. This procedure returns the absolute value of the number a. The result is
  1368. the variable a.
  1369.  
  1370. ****************************************************************************
  1371.  
  1372. PROCEDURE floorm(VAR a : mfloat);               { *** a <-- floor(a) *** }
  1373.  
  1374. The greatest integer number, which is less than or equal the mfloat number
  1375. a, is calculated and returned to the variable a.
  1376.  
  1377. ****************************************************************************
  1378.  
  1379. PROCEDURE fmodm(VAR a, b : mfloat);             { *** a <-- fmod(a,b) *** }
  1380.  
  1381. This is a floating point modul function. The number a is divided by b, so
  1382. that the result is an integer. The rest of this division is returned by the
  1383. variable a. The variable b is not changed.
  1384.  
  1385. a = k * b + c;   k is an integer;   0 <= c / b < 1;  c is assigned to a
  1386.  
  1387. If the variable b is zero, the error flag is set. For negative b, the result
  1388. is negativ or zero. For very large k, which can not be represented as an
  1389. mfloat number exactly, the result is zero.
  1390.  
  1391. ****************************************************************************
  1392.  
  1393. PROCEDURE frexpm(VAR a : mfloat;                { *** a <-- frexp(a,b) *** }
  1394.          VAR b : integer);
  1395.  
  1396. The mfloat number a is splitted in the exponent and in the mantissa:
  1397.  
  1398. a = c * (2**b);  0.5 <= |c| < 1 if a <> 0;  c is assigned to a
  1399.  
  1400. This procedure returns two results. IF the mfloat number is zero, it remains
  1401. zero and b is set to zero.
  1402.  
  1403. ****************************************************************************
  1404.  
  1405. PROCEDURE hypotm(VAR a, b : mfloat);            { *** a <-- hypot(a,b) *** }
  1406.  
  1407. This function is used for the calculation of the absolute value of a complex
  1408. number:
  1409.  
  1410. c = sqrt(sqr(a) + sqr(b)); c is assigned to a; b remains unchanged.
  1411. |a + i*b| = hypot(a,b);
  1412. This procedure calculates the correct result also if sqr(a) can not be
  1413. represented as a mfloat number, because it is too great or too small.
  1414.  
  1415. ****************************************************************************
  1416.  
  1417. PROCEDURE ldexpm(VAR a : mfloat; b : integer); { *** a <-- ldexp(a,b) *** }
  1418.  
  1419. The mfloat number a is multiplied by the power of two 2**b.
  1420.  
  1421. a <-- a * (2**b);
  1422.  
  1423. The calculation is very fast because only the exponent of two is manipulated.
  1424.  
  1425. ****************************************************************************
  1426.  
  1427. PROCEDURE logm(VAR a : mfloat);                   { *** a <- log(a) *** }
  1428.  
  1429. This procedure calculates the natural logarithm ln(x) of the argument a. The
  1430. result is the variable a.
  1431.  
  1432. a > 0                  range of the argument
  1433. -22712 <= a <= 22712   range of the result
  1434.  
  1435. ****************************************************************************
  1436.  
  1437. PROCEDURE log10m(VAR a : mfloat);              { *** a <- log10(a)  *** }
  1438.  
  1439. This procedure calculates the decade logarithm log(x) of the argument a. The
  1440. result is the variable a.
  1441.  
  1442. a > 0                  range of the argument
  1443. -9864 <= a <= 9863     range of the result
  1444.  
  1445. ****************************************************************************
  1446.  
  1447. PROCEDURE modfm(VAR a, b : mfloat);         { *** a, b <-- modf(a) *** }
  1448.  
  1449. This procedure splits the number a in an integer part and in an fractional
  1450. part:
  1451.  
  1452. a = k + b;    k is an integer;    0 <= b < 1;    k is assigned to a
  1453.  
  1454. ****************************************************************************
  1455.  
  1456. PROCEDURE powm(VAR a, b : mfloat);          { *** a <-- pow(a,b) *** }
  1457.  
  1458. This procedure calculates the power of the basis a by the exponent b.
  1459. The result is assigned to the variable a, the variable b is not changed.
  1460. If the exponent is an integer, the calculation is made by repeated
  1461. multiplications. The number of multiplications is proportional to
  1462. log(|b|). If b has a fractional, the calculation is made by logarithm. In
  1463. this case, the base a must be positiv, otherwise the error flag is set,
  1464. because the result is not real. If the number b is very large so that
  1465. no fractional can be represented by the mfloat format and the base
  1466. is negativ, the error flag is set too. If a = 0 and b <= 0 the result is
  1467. not defined and the error flag is set.
  1468.  
  1469. ****************************************************************************
  1470.  
  1471. PROCEDURE pow10m(VAR a : mfloat; b : integer);  { *** a <-- pow10(a) *** }
  1472.  
  1473. The power of the basis ten is build with the integer exponent b. The result
  1474. is assigned to a.
  1475.  
  1476. a <-- 10**b
  1477.  
  1478. if b < -9864 then a = 0
  1479. if b > 9863  then there is an overflow -> error flag is set
  1480.  
  1481. ****************************************************************************
  1482.  
  1483. PROCEDURE sinm(VAR a : mfloat);                    { *** a <- sin(a) *** }
  1484.  
  1485. This procedure calculates the sin(x) function of the argument a. The result
  1486. is the variable a.
  1487.  
  1488. ****************************************************************************
  1489.  
  1490. PROCEDURE sinhm(VAR a : mfloat);                  { *** a <- sinh(a) *** }
  1491.  
  1492. This procedure calculates the sinh(x) function of the argument a. The result
  1493. is the variable a.
  1494. -22712 < a < 22712   range of the argument without an overflow
  1495.  
  1496. ****************************************************************************
  1497.  
  1498. PROCEDURE sqrtm(VAR a : mfloat);                  { *** a <- sqrt(a) *** }
  1499.  
  1500. This procedure calculates the square root of the number a. The result is
  1501. variable a. If a is negativ, the error flag is set.
  1502.  
  1503. ****************************************************************************
  1504.  
  1505. PROCEDURE tanm(VAR a : mfloat);                   { *** a <- tan(a) *** }
  1506.  
  1507. This procedure calculates the tan(x) function of the argument a. The result
  1508. is the variable a. It is faster than the calculation of sin(a) / cos(a).
  1509.  
  1510. ****************************************************************************
  1511.  
  1512. PROCEDURE tanhm(VAR a : mfloat);                  { *** a <- tanh(a) *** }
  1513.  
  1514. This procedure calculates the tanh(x) function of the argument a. The result
  1515. is the variable a.
  1516.  
  1517. -1 <= a <= 1       range of the result
  1518.  
  1519. ****************************************************************************
  1520. ++++++++++++++++++    extended standart functions   ++++++++++++++++++++++++
  1521. ****************************************************************************
  1522.  
  1523. PROCEDURE acoshm(VAR a : mfloat);                { *** a <- arcosh(a) *** }
  1524.  
  1525. This procedure calculates the invers cosh(x) function of the argument a. The
  1526. result is the variable a.
  1527. 1 <= a          range of the argument
  1528. 0 <= a < 22712  range of the result
  1529.  
  1530. ****************************************************************************
  1531.  
  1532. PROCEDURE acotm(VAR a : mfloat);                 { *** a <- arccot(a) *** }
  1533.  
  1534. This procedure calculates the invers cot(x) function of the argument a. The
  1535. result is the variable a.
  1536. 0 < a < pi       range of the result
  1537.  
  1538. ****************************************************************************
  1539.  
  1540. PROCEDURE acothm(VAR a : mfloat);                { *** a <- arcoth(a) *** }
  1541.  
  1542. This procedure calculates the invers coth(x) function of the argument a. The
  1543. result is the variable a.
  1544.  -1 > a   or   1 < a    range of the argument
  1545.  
  1546. ****************************************************************************
  1547.  
  1548. PROCEDURE asinhm(VAR a : mfloat);                { *** a <- arsinh(a) *** }
  1549.  
  1550. This procedure calculates the invers sinh(x) function of the argument a. The
  1551. result is the variable a.
  1552.  
  1553. ****************************************************************************
  1554.  
  1555. PROCEDURE atanhm(VAR a : mfloat);                { *** a <- artanh(a) *** }
  1556.  
  1557. This procedure calculates the invers tanh(x) function of the argument a. The
  1558. result is the variable a.
  1559. -1 < a < 1  range of the argument
  1560.  
  1561. ****************************************************************************
  1562.  
  1563. PROCEDURE cossinm(VAR a, b : mfloat); { *** a <-- cos(a), b <-- sin(a) *** }
  1564.  
  1565. This procedure calculates the sin(x) and cos(x) function by only one
  1566. evulation of the series for the sine function. This procedure is therefore
  1567. faster than two individual calculations. It is recommended for example
  1568. for the calculation of the exponential function of complex arguments.
  1569.  
  1570. ****************************************************************************
  1571.  
  1572. PROCEDURE cotm(VAR a : mfloat);                     { *** a <- cot(a) *** }
  1573.  
  1574. This procedure calculates the cot(x) function of the argument a. The result
  1575. is the variable a. It is faster than the calculation of cos(a) / sin(a).
  1576.  
  1577. ****************************************************************************
  1578.  
  1579. PROCEDURE cothm(VAR a : mfloat);                   { *** a <- coth(a) *** }
  1580.  
  1581. This procedure calculates the coth(x) function of the argument a. The result
  1582. is the variable a.
  1583. range of the result: a >=1 or a <= -1
  1584.  
  1585. ****************************************************************************
  1586.  
  1587. PROCEDURE exp10m(VAR a : mfloat);                  { *** a <- 10 ** a *** }
  1588.  
  1589. This procedure calculates the power of 10 of the argument a. The
  1590. result is the variable a.
  1591.  a < 9863    range of the argument without an overflow
  1592.  a < -9864   there is an underflow, result = 0
  1593.  a >= 0      range of the result
  1594.  
  1595. ****************************************************************************
  1596.  
  1597. PROCEDURE sqrm(VAR a : mfloat);                     { *** a <- sqr(a) *** }
  1598.  
  1599. This procedure calculates the sqare of a. The result is variabel a.
  1600. sqr(a) = a * a
  1601.  
  1602. ****************************************************************************
  1603.  
  1604. PROCEDURE truncm(VAR a : mfloat);              { *** a <-- trunc(a) *** }
  1605.  
  1606. This procedure returns the integer part of the number a (truncation). The
  1607. result is variable a. The sign has no influence and remains unchanged.
  1608.  
  1609.  
  1610. ===========================================================================
  1611. ***************************************************************************
  1612. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  1613. ***************************************************************************
  1614. ===========================================================================
  1615.  
  1616. Examples with MFLOAT:
  1617.  
  1618. a) The calculation of the number PI:
  1619.    There are many methods known to calculate the irrational number pi,
  1620.    which is the relation between the perimeter and the diameter of a circle.
  1621.  
  1622.    Pi = 16 * arctan(1 / 5) - 4 * arctan(1 / 239);
  1623.    Pi = 48 * arctan(1 / 18)  + 32 * arctan(1 / 57) - 20 * arctan(1 / 239);
  1624.    Pi = 6  * arctan(1 / sqrt(3));
  1625.    PI = 4 * (1 - 1/3 + 1/5 - 1/7 + ....)
  1626.  
  1627.    The function arctan(x) is calculated with the series:
  1628.  
  1629.    arctan(x) = x - x**3 / 3 + x**5 / 5 - x**7 / 7 + ......
  1630.  
  1631.    The operation "x**.." means the power of x.
  1632.  
  1633.    The advantage of the first two equations is, that they converge very
  1634.    fast and every element of the series can be calculated recursivly from
  1635.    the former element by only two divisions by an integer. A division by
  1636.    an integer is very fast for long MFLOAT-numbers.
  1637.    The third equation is build by the known value : tan(pi/6) = 1 / sqrt(3).
  1638.    This method is slower than the former. The forth equation has only
  1639.    theoretical importance, in the numerical use this equation is
  1640.    very bad.
  1641.    The example program compares the different methods to calculate pi.
  1642.  
  1643. b) The evulation of the bessel function for large arguments:
  1644.    There are some series known, where some elements are very large and
  1645.    the result is small. In this case numerical errors are very large.
  1646.    An example is the bessel function J0(x) for an argument x = 100.0.
  1647.  
  1648.    J0(x) = 1 - (x/2)**2 / (1! * 1!) + (x/2)**4 / (2! * 2!) - ......
  1649.  
  1650.    The series converges very good, but there are problems with the accuracy
  1651.    of the calculation. The test program shows the influence of the accuracy
  1652.    of the calculation to the result.
  1653.  
  1654. c) Conversion of cartesian coordinates to polar coordinates:
  1655.    This example shows the use of the functions hypot(x,y) and atan2(x,y).
  1656.    This program also demonstrates the input of MFLOAT-numbers.
  1657.  
  1658. ===========================================================================
  1659. ***************************************************************************
  1660. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  1661. ***************************************************************************
  1662. ===========================================================================
  1663.  
  1664. Graz, 28.4.1993
  1665.  
  1666. Kaufmann Friedrich & Mueller Walter
  1667. Students at the Technical University of Graz
  1668.  
  1669. Registration address:
  1670.  
  1671. Kaufmann Friedrich
  1672. Schuetzenhofgasse 22
  1673. A-8010 GRAZ
  1674. AUSTRIA
  1675. EUROPE
  1676.  
  1677. email address:
  1678. fkauf@fstgds06.tu-graz.ac.at
  1679.